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1 Introduction 

Markov Chain Monte Carlo methods are statistical 
tools that have been recently proposed for the 
resolution of the MEG inverse problem [1]. Their main 
advantages are easy incorporation of a priori 
knowledge, and an adequate response to the ambiguity 
of the ill-posed MEG inverse problem. 

However, since simpler MCMC schemes might have 
difficulties in escaping from local modes, adequate 
research is needed to assess the most efficient MCMC 
algorithms for the resolution of the MEG inverse 
problem [2,3]. We present here both our results 
concerning MCMC scheme comparison, and the 
imaging tools we have developed for source 
computation and display. 

We first describe the probabilistic MEG model and the 
MCMC algorithms we have used (Parallel Tempering 
[4] and Reversible Jump [5]). In a first simulation 
study with a spherical head model, the combination of 
these MCMC schemes is shown to be superior to 
Simulated Annealing, both in accuracy and in its 
ability to estimate an unknown number of sources. 

We then apply Parallel Tempering to simulated data in 
a real head model obtained from MRI data. We 
describe how the MRI data are automatically processed 
to extract the brain, so that the source positions can be 
restricted to the cortex layer. We then compare Parallel 
Tempering with Simulated Annealing and show that 
PT exhibits a better robustness to noise level. 

We then use the 3-dimensional information previously 
extracted from the MRI images to display the MCMC 
results with this 3-dimensional morphological 
information. This further demonstrates how MCMC 
answers the MEG inverse problem by giving it a 
probabilistic solution. 

2 Methods 

2.1 Bayesian model and MCMC schemes 

Using a Bayesian model, the posterior probability of 
the sources is given by the product of the likelihood 
and the prior probability: 


P(J|B*)«p(B*|J).p(J) (1) 


where B 0 b S is the set of the M observed measurements 


(bi, b 2 , b M ), and J = (ji, j 2 , Jn) is a vector of N 
source dipoles. The prior probability p(J) enables easy 
incorporation of a priori information on the sources 
position, strength, and also number if this parameter is 
a priori unknown. Our analysis is limited to a single 
time slice, so that there is no dependence of the data 
and parameters on time. The likelihood is computed 
under the assumption of Gaussian noise. 

Finding a single value that maximizes the posterior 
only partly solves the MEG inverse problem, since it 
doesn’t convey any information on other possible 
sources. Therefore, we estimated the distribution of the 
posterior probability by sampling it with MCMC 
algorithms : we used Parallel Tempering [4] to prevent 
convergence to local modes, and Reversible Jump [5] 
to estimate the number of sources when this parameter 
was considered unknown. 

Parallel Tempering (PT) consists in running several 
Markov chains Q with probability distributions 
defined by : 

/! = (p(J|B*)Pi = l-k (2) 

To compute the (n+l) th realizations, two moves are 
possible : one can either conduct k independent 

Metropolis moves, or swap the values J- n) and J-^ of 
2 adjacent chains with the acceptance ratio : 
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Chains with high temperature ensure a fast and large 
exploration of the solution space, while chains with 
smaller temperature focus on the high probability 
regions. 

When the number of sources was a priori unknown, we 
combined PT algorithm with Reversible Jump. In this 
case, the solutions take the form (n, J n ), where n is the 
estimated number of sources, and J n is a n-source 


vector. The solution space is extended to |^J J (l) x {i}, 

i 

where J (l) is the set of all the i-source vectors. 



Reversible jumps permit movement between solutions 
of different dimensions by adding/deleting a random 
source, randomly splitting a source into 2 dipoles, or 
merging 2 sources into 1. These moves from dimension 
i to j are accepted with probability : 
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This results in a final algorithm that combines 3 
different moves : independent Metropolis update, PT 
moves and reversible jumps. 


2.2 MRI data processing 

For real data and simulations with a realistic head 
model, T1-weighted MRI axial images were 
acquired (voxel size = 0.9375 2 x 1.5 mm 3 ), and 
automatically processed as in [6]. First, the total 
histogram of the MRI data is smoothed with 
Silverman bump-hunting method [7], in which the 
grey level histogram is estimated by 
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(5) 


where N is the number of voxels, f is a Gaussian 
function N(0,1) and h is a parameter controlling the 
degree of smoothing. Beginning with large values of 
h, our algorithm makes a dichotomic search for the 
smallest value that will produce a distribution with 
exactly 4 local maxima (background, head, grey and 
white matter). The smoothed histogram is then 
automatically thresholded by searching points with a 
null derivative. 

Since there is no unambiguous correspondence 
between MRI grey levels and anatomical structures, 
this simple segmentation must still be refined. Using 
mathematical morphology tools such as erosion and 
dilation [8], we first separate the whole head object 
from the background. The brain is then correctly 
extracted by a series of erosions to remove spurious 
brain areas, and dilations to reconstruct the whole 
brain object (see Figure 1). For 3D visualization, the 
boundaries of the brain and head are computed with 
the Marching Cube algorithm [9] modified to allow 
for variable level of details, and to resolve the 
ambiguities of the initial algorithm [10]. 


3 Results 

3.1 Spherical model simulation 

A first simulation study was conducted to compare 
the combination of Parallel Tempering and 
Reversible Jump to Simulated Annealing. The head 
was modeled by a sphere of radius 80 mm. 


Measurements for 2 dipoles were simulated for 80 
radial magnetometers, with both white noise (20dB 
SNR) and neuromagnetic noise (60 random sources, 
20 dB SNR). For the prior, we used a 1 to 100 ratio 
for the inner/outer 40 mm shells of the brain sphere. 
The angle between the position vector and the 
magnetic moment followed a normal law N(7e/2, 
7t/6), and the number of sources was modeled by a 
Poisson prior between 1 and 4. Convergence of the 
MCMC sampler was monitored with the posterior 
log-probability and the percentage of explained 
variance (see Figure 2). 



Figure 1: Result of the MRI segmentation and 
mathematical morphological filtering, showing the 
extracted brain and head objects. This treatment 
was performed automatically, in less than 30 
seconds on a standard PC computer (124 slices of 
256 x 256 voxels). 



Figure 2: Explained variance, against the iteration 
number. Convergence of the chain is obtained when 
explained variance and posterior probability 
become stable. The PT/RJ algorithm was run for 
3000 iterations ; the average number of iterations 
required for convergence was less than 1400. 


Table 1: Mean error for source position, showing the 
better accuracy obtained with PT/RJ, even with the 
added difficulty of an unknown number of dipoles. 



Iterations 

Sources 

Mean error 

PT/RJ 

3.000 

1-4 

4.5 mm 

SA 

10.000 

2 

6.3 mm 





















30 trials were conducted for the PT/RJ algorithm and 
for Simulated Annealing. It was found that the mean 
error for source position was smaller for PT/RJ than for 
SA (see Table 1). Moreover, whereas the number of 
sources must be a priori known for SA, the PT/RJ 
algorithm was able to reliably estimate the true number 
of sources (see Figure 3). 
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Figure 3 : Estimated number of dipoles, showing a high 
probability for a 2-source pattern. 

3.2 Realistic model simulation 

A second simulation was conducted to study the 
influence of noise on the Parallel Tempering 
algorithm. MRI images from one subject were 
processed as described in 2-2, and measurements for 
2 dipoles were simulated with both white noise 
(20dB SNR) and neuromagnetic noise (5, 10, 15 and 
20 dB SNR). Sources were restricted to the brain 
volume, and constrained to the cortex (approximated 
by the 5 outer mm of the brain) by a 1 to 1000 
probability ratio. Other priors were as in the 
previous study. 

For each noise condition, 20 runs were made for PT 
(30.000 iterations) and for SA (150.000 iterations). 
Only the runs that explained at least 95% of the 
measurement variance were considered as valid and 
were used for computing the mean error in source 
position. The other runs were rejected without 
further analysis. 

Although SA and PT gave equivalent results for low 
noise level (16 dB SNR), SA results degraded 
quickly with increasing noise level (45% of valid 
runs for 8 dB SNR). On the contrary, PT was little 
affected by noise level, and it was only for very high 
noise level (4.1 dB SNR) that PT performance 
degraded very quickly (see Table 2). And we feel 
that in this case, the failure of PT was more due to 
the incapacity of explaining the measurements with 
our 2-source model, than from erroneous 
convergence to local modes. 

To display MCMC results, source position were 
discretized according to the MRI voxel grid, and for 
each MRI voxel, we computed and displayed the mean 
intensity of MCMC solutions. The results clearly 


showed the position of the 2 source solutions, as well 
as their likely extent (see Figure 4). 

Table 2 : influence of noise on PT and SA results. 
Although PT and SA give similar results for low 
noise, PT was much more robust to increasing noise 
level than SA. 


SNR 

SA 

%valid 

PT 

%valid 

SA 

error 

PT 

error 

15.8 dB 

100% 

90% 

8.7 

7.0 

10.1 dB 

70% 

90% 

9.5 

9.4 

7.8 dB 

45% 

95% 

10.6 

11.5 

4.1 dB 

0% 

15% 

- 

17.7 




: Dipole distribution obtained with PT 
algorithm (SNR = 8 dB). Note the information 
provided by MCMC for both source position and 
likely extent. 


4 Discussion 

Our results show that adequately chosen MCMC 
algorithms are able to solve the MEG inverse 
problem even under difficult conditions like an 
unknown number of sources (RJ + PT algorithm), or 
a high noise level (PT). Moreover, MCMC methods 
have proven to be superior to SA in terms of 






























accuracy, estimation of source number, and 
robustness to noise. In addition, MCMC methods 
give an estimate of the source probability 
distribution, so that we not only have a single best 
solution, but a set of likely solutions, as well as their 
extent. This last point has already been stressed in 
[1], and we feel that this may constitute one of the 
most appealing advantages of MCMC methods. 
However, a number of issues have still to be solved 
for practical application of MCMC methods. The 
most important one is undoubtedly their extension to 
spatio-temporal analysis, a point that is already 
covered by others in [11]. In addition, MCMC 
results by themselves provide little or no 
information on whether they encompass the whole 
set of possible solutions, or only a part of it. This 
point has already been studied in a number of both 
theoretical and practical papers on MCMC 
convergence assessment, and undoubtedly, the 
reliability of MEG inverse problem resolution would 
greatly benefit from the application of such tools. 
Last, MCMC pose great challenges for accurate 
graphic display of their solutions, so that research is 
also needed for an efficient simplification and 
display of the posterior probability estimate. 
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